function DecodingProbability_v16()

% Clear workspace
clear all
close all
clc

% Time it
tic

% Transmitter range
T_min=20;
T_step=1;
T_max=60;

% Field cardinality
q=2^8;

% Layer selection probability
g=[0.3 0.7]; % 0.5 0.5

% Layer dimensions
k1=15; % 33
k2=25; % 100-33
K1=k1;
K2=k1+k2;

% Decoding 1. Layer probabilities
l1_prob=zeros(T_max,1);

for pkts_recv = T_min:T_step:T_max % For each nb of transmitted packets
    l1_prob(pkts_recv)=MatrixFun(pkts_recv,k1,k1,g(1),q);
    disp(['Layer 1: ' num2str(pkts_recv) ' out of ' num2str(T_max)])
end

% Decoding 2. Layer probabilities
l2_prob=zeros(T_max,1);

% Array for temporary results
sol=zeros(T_max+1,1); 

for tx =T_min:T_step:T_max % For each number of recv packets
    
    nn=1; % matlab indexing workaround
    
    binom_vector=binopdf(0:tx,tx,g(1));
    
    for n = 0:tx % For all permutations of recv packets
        
        val=0;
        
        for i=0:K1 % For all possible ways to achieve rank K2 with given permutation of recv packets
            tmp1=ProbMatricesWithRank(n,K1,i,q);
            tmp2=ProbMatricesWithRank(tx-n,K2-i,K2-i,q);
            val=val+tmp1*tmp2;
        end
        
        sol(nn)=val*binom_vector(nn);
        nn=nn+1;
        
    end
    
    l2_prob(tx)=sum(sol);
    sol=zeros(T_max,1);
    
    disp(['Layer 2: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

% answer








% Wait, we know that whenever we have layer 2, we also have layer 1
%
%


% % new_l1_prob=zeros(T_max,1);
% % for i =K2:length(l1_prob)
% % %     new_l1_prob(i)=2-l2_prob(i)-(1-l2_prob(i))/(1-l1_prob(i)); % try 1
% % %     new_l1_prob(i)=l2_prob(i)/(l2_prob(i)+(1-l2_prob(i))*l1_prob(i)); % try 2
% %     new_l1_prob(i)=l2_prob(i)*l1_prob(i)/(l2_prob(i)+l1_prob(i)); % try 2
% % end
% % 
% % % l1_prob
% % % new_l1_prob
% % 
% % l1_prob(K2:T_max)=new_l1_prob(K2:T_max);
% % 
% % 
% % 




plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2)

% Time it
toc

end

% Plotter for a nice graph!
function plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2)

% Replace 0 with NaN in (l1_prob,l2_prob) for prettier plot
for k=1:length(l1_prob)
    
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
    
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
    
end

% Plotting
figure(1)
hold('on')
plot(1:K2-1,l1_prob(1:K2-1,1),'-*','Color','b')
plot(1:T_max,l2_prob,'-*','Color','g')
hold('off')

% Plot annotation
legend('Layer 1','Layer 2','location','SouthEast')
grid('on')
% pbaspect([2.5 1 1])
set(gca,'XTick',0:10:T_max)
set(gca,'YTick',0:0.1:1)
xlim([T_min T_max])
ylim([0 1])

% Save plot
% print(gcf,'uep_ew_analytic.eps')

end

% Shady function, eliminate it!
function Pr = MatrixFun(pkts_recv,layer_length,rank,g,q)
% This function returns the probability of a matrix reaching rank 'rank'
% for a number of 'pkts_recv', where there is a prob. of 'g' it being to
% for this layer.

Pr=0;

if pkts_recv>=rank % if we have less pkts_recv than required rank, no need to calculate -> it is = 0
    
    bino_vector=binopdf(1:pkts_recv,pkts_recv,g); % binomial vector with the probability of 1,2,3,...,pkts_recv from layer with prob 'g'
    
    for outcome=1:length(bino_vector) % for each outcome of the received packets
        if outcome>=rank
            Pr=Pr+bino_vector(outcome)*ProbMatricesWithRank(outcome,layer_length,rank,q);
            %             [bino_vector(outcome) outcome layer_length rank]
        end
    end
    
else
    Pr=0;
end

assert(isnan(Pr)==0,['MatrixFun returned NaN with:' num2str([pkts_recv,layer_length,rank,g,q])]);

end

% Prob. of matrices with dimensions (m,n) of rank (r) over field size (q)
function PMWR = ProbMatricesWithRank(m,n,r,q)

% The probability og getting rank 0 is = 0 for m=r=0
if m==0 && r==0
    PMWR=1;
    return
end

if n==0
    'RØØØØØØV'
end

% If m<r there is no possibility to achieve r rank
if m<r
    PMWR=0;
    return
end

% Get first set of gaussian coefficients
gc=gausscoeffs2(n,r,q);

% % % Calculate "sum"
% % val=0;
% % for k=0:r
% %     % This should be the one!
% %     val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(m*k+binomcoeffs(r-k,2)-n*m));
% % end

% Calculate "sum"
val=0;
for k=0:r
    % This should be the one!
    val1=gausscoeffs2(r,k,q);
    val2=q^(m*k+binomcoeffs(r-k,2)-n*m);
    
    %     assert(isnan(val1)==0,'val1')
    %     assert(isnan(val2)==0,'val2')
    
    % Multiply by 0 before inf :)
    if val1==0 || val2==0
        comb_val=0;
    else
        comb_val=val1*val2;
    end
    
    val=val+((-1)^(r-k)*comb_val);
    
    %     assert(isnan(val)==0,'ProbMatricesWithRank returned NaN');
    
end

% Return probability of matrix 'm'x'n' with rank 'r'
PMWR=gc*val;

assert(isnan(PMWR)==0,['ProbMatricesWithRank returned NaN, params: ' num2str([m n r q])]);

end

% Compute gaussian coefficients (See Wolfram Alpha)
% Version 2: Numerical improvements over gausscoeffs(... )
function GC = gausscoeffs2(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    % There is X factors in numerator
    
    
    num_length=abs((m))-abs((m-r+1));
    num=ones(num_length,1);
    num_index=1;
    
    for w=m:-1:m-r+1
        num(num_index,1)=(q^w-1);
        num_index=num_index+1;
    end
    
    % Calculate denominator
    %     denom=1;
    
    denom_length=r;
    denom=ones(denom_length,1);
    denom_index=1;
    
    for w=r:-1:1
        denom(denom_index,1)=(q^w-1);
        denom_index=denom_index+1;
    end
    
    % Calculate gaussian coefficient
    GC=prod(num(:)./denom(:));
    
    assert(isnan(GC)==0,'Exception: gausscoeffs2 returned NaN')
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end

end

% Binomial polynomial thing
% As on page 123 in "A course in combinatorics"
function BC = binomcoeffs(a,k)

tmp_vector=ones(2,1);
tmp_index=1;

for w=0:-1:-k+1
    tmp_vector(tmp_index)=(a+w);
    tmp_index=tmp_index+1;
end

num=prod(tmp_vector);
denom=factorial(k);

BC=num/denom;

end





















%% Deprecated

% Should work (Tested! see bottom)
% Deprecated

% function GC = gausscoeffs(m,r,q)
% if r==0
%     % disp('r = 0 in gauss coeffs')
%     GC=1;
% elseif r>0
%     % disp('r > 0 in gauss coeffs')
%
%     % Calculate numerator
%     num=1;
%     for w=m:-1:m-r+1
%         num=num*(q^w-1);
%     end
%
%     % Calculate denominator
%     denom=1;
%     for w=r:-1:1
%         denom=denom*(q^w-1);
%     end
%
%     % Calculate gaussian coefficient
%     GC=num/denom;
%
% elseif r<0
%     disp('r < 0 error in gausscoeffs!!!')
% end
%
%
% end


